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We investigate short range order in liquid and supercooled liquid Fe and Fe-based metallic glass 
using ab-initio simulation methods. We analyze the data to quantify the degree of local icosahedral 
and polytetrahedral order and to understand the role of alloying in controlling the degree of geometric 
frustration. Comparing elemental Fe to Cu l| we find that the degree of icosahedral order is greater 
in Fe than in Cu, possibly because icosahedral disclination line defects are more easily incorporated 
into BCC environments than FCC. In Fe-based metallic glass- forming alloys (FeB and FeZrB) we find 
that introducing small concentrations of small B atoms and large Zr atoms controls the frustration 
of local icosahedral order. 

PACS numbers: 61.43.Dq,61.20. Ja,61.25.Mv 

I. INTRODUCTION 

As noted by Frank 2|, the local icosahedral clustering of 12 atoms about a sphere is energetically preferred because 
it is made up entirely of four-atom tetrahedra, the densest-packed cluster possible. However, local icosahedral order 
cannot be propagated throughout space without introducing defects. Frustration of packing icosahedra is relieved in 



a curved space, where a perfect 12-coordinated icosahedral packing exists 3, |4|, l5|. We adopt this structure as an 



ideal reference against which actual configurations will be compared. 

Disclination line defects of type ±72° may be introduced into this icosahedral crystal and thereby control the 
curvature. In order to "flatten" the structure and embed it in ordinary three dimensional space an excess of —72° 
disclinations is needed, and these cause increased coordination numbers of 14, 15 or 16. Large atoms, if present, 
would naturally assume high coordination numbers and aid in the formation of a disclination line network. Similarly, 
smaller atoms would naturally assume low coordination numbers of 8, 9 or 10, and have positive disclinations attached 
to them, increasing the frustration. For a particular coordination number, it may be possible to construct a cluster, 



known as a Kasper-Polyhedron (see Section UlI C|) . made entirely of tetrahedrons. 



Honeycutt and Andersen 



6[ introduced a method to count the number of tetrahedra surrounding an interatomic 



bond. This number is 5 for icosahedral order with no disclination, 6 for a —72° disclination and 4 for a +72° 
disclination. Steinhardt, Nelson and Ronchetti 7[ introduced the orientational order parameter Wq to demonstrate 



short range icosahedral order. We employ both methods to analyze icosahedral order in supercooled metals and metal 



alloys, in addition to conventional radial distribution functions, structure factors and Voronoi analysis. 



Many simulations have been performed on pure elemental metals and metal alloys using model potentials, 8l.l9l, 1 10 



11], but do not necessarily produce reliable structures owing to their imperfect description of interatomic interactions. 



First principles {ab-initio) calculations achieve the most realistic possible structures, unhindered by the intrinsic 
inaccuracy of phenomenological potentials and with the ability to accurately capture the chemical natures of different 
elements and alloys. The trade-off for increased accuracy is a decrease in the system sizes one can study, so only local 



order can be observed, not long range. Also runs are limited to short time scales. Early ab-initio studies on liquid 



studies on Ni and Zr 
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Copper [12], [13, |14[ and Iron [15j were not analyzed from the perspective of icosahedral ordering. Recent ab-initio 



17j find that the degree of icosahedral ordering increases with supercooling in Ni, while in 



Zr BCC is more favored. Studies on binary metal alloys by Jakse et al. 18|, ll9( and by Sheng et al 20| quantify local 



icosahedral order in the alloys. We previously [lj investigated icosahedral order in liquid and supercooled Cu. 

Elemental metals crystallize so easily that they can hardly be made amorphous at any quench rate. Alloying can 
improve the ease of glass formation. For some special alloys, a bulk amorphous state can be reached by slow cooling. 
Pure elemental Fe is a poor glass former, but Fe-based compounds like FeB and especially FeZrB, show improved glass 



formability. We augment our molecular-dynamics simulation with another algorithm called 'Tempering' or 'replica 



exchange method' (REM) 2jJ, l22j for fast equilibration at low temperatures 



In comparison to liquid and supercooled liquid Copper [1| which show only weak icosahedral order and very little 



temperature variation, Fe showed a monotonic increase in icosahedral order, which became very pronounced when 
supercooled. Analysis of quenched Fe revealed a natural way of introducing 5-fold coordinated bonds plus a single -72° 
disclination line segment into an otherwise perfect BCC environment, without disturbing the surrounding structure. 
Addition of B to Fe decreased the icosahedral order, due to the positive disclinations centered on the smaller B atom 



which increased frustration. Further inclusion of larger Zr atoms to form FeZrB found an enhanced icosahedral order 



compared to FeB. This could possibly be explained by formation of negative disclination line defects 23| anchored on 



the larger Zr atoms, which eases the frustration of icosahedral order on the Fe atoms. 



At high temperatures all of our measured structural properties of liquid Cu |l[ and liquid Fe resembled each other, 



and also resembled a maximally random jammed [2J] hard sphere configuration. This suggests that a nearly universal 



structure exists for systems whose energetics are dominated by repulsive central forces. 
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Section |TT] describes our combined method of monte-carlo and first principles MD, that we refer to as "Tempering 
MD" and discusses other simulation details. Section HTT1 presents our results on pure Fe while fsection HV|) compares 
this with FeB and FeZrB alloys. 

II. TEMPERING MOLECULAR DYNAMICS (TMD) AND OTHER SIMULATION DETAILS 

One reason alloys for m g lass more easily is that chemical identity introduces a new configurational degree of freedom 



that evolves slowly 
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261 ] . Unfortunately, this makes simulation more difficult. It is especially difficult to equilibrate 



the system at very low temperatures, because the probability to cross an energy barrier drops, trapping it in particular 



configurations. For this reason we use a Monte-Carlo method, known as tempering or replica exchange 2JJ, [22| to 



augment our first-principles MD, allowing us to sample the configurational space more efficiently than conventional 
MD. 

In the canonical ensemble, energy fluctuates at fixed temperature. A given configuration C with energy E can occur 
at any temperature T with probability proportional to e~ l3E , (j3 — 1/kT). Now consider a pair of configurations, C\ 
and C'2 of energy E\ and E2 occurring in simulations at temperatures T\ and 2V We can take C\ as a member of the 
ensemble at T2 , and C2 as a member of the ensemble at T\ , with a probability 

p = e -(P2-Pi)(Ei-E 2 ) /js 

without disturbing the temperature-dependent probability distributions of energy (or any other equilibrium property) . 
Because each run remains in equilibrium at all times even though its temperature changes, we effectively simulate a 
vanishingly low quench rate. 
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TABLE I: Details of tempering MD runs. 



Chemical Species 


Temperatures(K) 


density (A -3 ) 


time (ps) 


Feioo 


800-1900 


0.0756 


1.5 


Fe8oB2o 


800-1500 


0.0814 


1.8 


Fe7oZrioB2o 


800-1800 


0.0787 


1.8 



In practice we perform several MD simulations at temperatures separated by 100K. We use ultrasoft pseudopoten- 



tials [27| as provided with VASP 28j to perform the MD simulation. All calculations are T' point calculations (a 
single 'k' point). All runs use an MD time step of 2 fs, and reach total simulated time of order 1.5-1.8 ps (see Table|J) 
with a total of N=100 atoms. Every 10 MD steps we compare the energies of configurations at adjacent temperatures 
and swap them with the above probability. Eventually, configurations initially frozen at low temperature reach a 
higher temperature. The simulations then can carry the structure over energy barriers, after which the temperature 
can again drop. 

In an effort to explore the structures of compounds with differing glass-forming ability we compare pure elemental 
Iron and two Iron-based glass- forming alloys. Tempering MD requires that we perform simulations at a constant 
density for all the temperatures, but we have no rigorous means of predicting the density at high temperature. For 



pure liquid Iron, the density is known experimentally [29j, and we use this value. For FeB and FeZrB, we took a 



high temperature liquid structure and quenched it, relaxing positions and cell lattice parameters, to predict a low 
temperature density. We then decreased the density of the relaxed structure by 6 percent to account for volume 
expansion, to arrive at the densities used in our liquid simulations. 

Because of the efficient sampling of our tempering MD method, the structure of pure Fe partially crystallizes at low 
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FIG. 1: Radial distribution function of pure elemental liquid Fe. Simulations (magnetic and non-magnetic) are run at T=1800K 



compared with X-ray experiments at T=1833K 



29( and Neutron scattering experiments at 1830K 



temperatures after about 1 ps. In the following discussion of our T=800K sample we will refer to different structural 
features before and after crystallization. We also performed several long (2. Ops) conventional first-principles MD at 
T=800K yielding results similar to the results of tempering MD prior to crystallization. 

For all runs we employed spin polarization, reasoning that local magnetic moments exist even above the Curie point. 
These local moments have a significant influence on the short-range order because ferromagnetic Iron prefers a longer 
bond length than paramagnetic Iron [30j] . Of course, the ferromagnetic state of the liquid implies improper long-range 
correlations. Unfortunately, since our forces are calculated for electronic ground states, we cannot rigorously model 



the true paramagnetic state of liquid Iron and Iron-based alloys with these methods. 



III. PURE FE 
A. Radial Distribution Function g(r) 

The radial distribution function, g(r), is proportional to the density of atoms at a distance r from another atom 
and is calculated here by forming a histogram of bond lengths. We use the repeated image method to obtain the 
bond lengths greater than half the box size and anticipate g(r) in this range may be influenced by finite size effects. 
Further, we smooth out the histogram with a gaussian of standard deviation 0.05A. 

To evaluate the role of magnetism, Fig. [1] illustrates the radial distribution functions for liquid Fe simulated at 
T=1800K, just below melting (T m =1833K). Evidently, the simulation with magnetism yields good agreement in the 
position and height of the first peak in g(r) with experimental X-ray g(r) [29(, while neglect of magnetic moments 
results in near neighbor bonds that are too short and too weak. However, magnetism overestimates the strength 
of long-range correlations beyond the nearest-neighbor peak, while neglecting magnetism yields reasonably accurate 
long-range g{r). Nevertheless, for the present study of local order, it is necessary to make spin polarized calculations, 
to get the short range correlations and hence the local order correct. Strangely, a recent experimental neutron 



g(r) 3l| has a shorter and broader first peak compared to both our magnetic g{r) and to the g(r) from the prior 



X-ray diffraction experiment 29|. The positions of the different maxima and minima in our simulated magnetic g(r) 
compare well with both the experiments. The position of the first peak in our magnetic g(r) is shifted by 0.05A to 
the left of the neutron first peak. The X-ray experiment doesn't have enough data points around the first maximum 
to determine the peak position accurately. 

We calculate the coordination number by counting the number of atoms within a cutoff distance from a central 
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atom. We choose the cutoff distance (R cu t) at the first minimum of g(r). For pure Fe the minimum is at R cu t=3.5 A. 
The precise location of the minimum is difficult to locate, and its variation with temperature is smaller than the error 
in locating its position, so that we don't change the value of R cu t with temperature. With this value of R cu t we find 
an average coordination number (N c ) of 13.2 which is nearly independent of temperature (N c changes from 13.1 at 
high temperature to 13.3 with supercooling). 

B. Liquid Structure Factor S(q) 

The liquid structure factor S(q) is related to the radial distribution function g(r) of a liquid with density p by, 

oo 

S(q) = 1 + Airp [ \g(r) - 1] ^Mr 2 dr (2) 

J qr 

o 

One needs the radial distribution function up to large values of r to get a good S(q). In our first principles simulation, 
we are restricted to small values of r, due to our small system sizes, so we need a method to get S(q) from our limited 



g(r) function. Baxter developed a method 32j, l33j to extend g(r) beyond the size of the simulation cell. The method 



exploits the short range nature of the direct correlation function c(r) , which has a range similar to the interatomic 



interactions [3J|, as opposed to g(r) which is long ranged. 

Assuming that c(r) vanishes beyond a certain cutoff distance r c , we solve the Baxter's equations iteratively to obtain 
the full direct correlation function for < r < r c . From c(r) we calculate the structure factor S(q) by a standard 
Fourier Transform. The S(q) showed good convergence with different choices of r c . A choice of r c =5A seemed 
appropriate because it was one half of our smallest simulation cell edge length. Even though in metals there are long 
range oscillatory Friedel oscillations, our ability to truncate c(r) at r c =5A shows that these are weak compared with 



short range interactions. An application of this method to obtain S(q) of Cu 1[ showed excellent agreement with the 



experimental S(q). 



The simulated S(q) for pure Fe at T=1800K (see Fig. ^ is compared to recent neutron scattering experiments at 
T=1830K. Even though the positions of the different peaks compare very well, there is serious discrepancy in their 
heights. Especially, the height of the first peak of our simulated S(q) is higher than that of the experiments. This 
discrepancy is expected because we include magnetism, which gives accurate short-range correlations while overesti- 



mating the long-range ones (see Fig [IT) . But the cause of the discrepency is not entirely clear since a comparison of 



the simulated structure factor of Ni 16[ (done without including magnetism) with neutron scattering experiments 31 1 



shows similar discrepancies between their SjqYs 



A sum rule can be obtained for S (q) 



35 



36|. By inverting the fourier transform of Eq.[2]and then taking the r — > 



limit, one gets 



Q 

I(Q) = J q 2 [S(q)-l}dq^-27T 2 p 



(3) 



in the limit Q — > oo. Further, the integral is supposed to oscillate with Q about the limiting value as Q — > oo. Using 



our S(q) we observed that the integral is consistent with the sum rule and oscillates nicely about the limiting value for 



Q > 3 A- 1 , while using the S(q) from the neutron scattering experiments 



3 11 ] , we observe a positive drift in the mean 



value about which the integral oscillates. Such a drift could indicate the presence of spurious background corrections. 



The S(q) from the X-ray experiment 29( seems to be in good agreement with the ideal sum rule. 



As we lower the temperature, the peak heights in S(q) grow, indicating an increase in short range order with 



supercooling. We also observe a slight shoulder in the second peak of the S(q) (Fig. [2]), which grows with supercooling. 
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FIG. 2: Comparison of simulated and experimental S(q) near the melting temperature of Fe. 
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FIG. 3: Distribution of We in liquid and supercooled liquid Fe and Cu. Fe shows more pronounced icosahedral order than Cu 



with supercooling 



The split second peak positions are in the ratios of 20:12 and 24:12 with respect to the first peak positions, just what 



one would ideally observe if there was icosahedral order 



23, 
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TABLE II: We values for a few clusters. The Z8-Z16 series are Kasper-Polyhedrons including the perfect icosahedron (Z12). 



Cluster 


No. of atoms 


Voronoi type 


We 


HCP 


12 


(0,12,0) 


-0.012 


FCC 


12 


(0,12,0) 


-0.013 


BCC 


14 


(0,6,0,8) 


+0.013 


Z8 


8 


(0,4,4) 


+0.010 


Z9 


9 


(0,3,6) 


-0.038 


Z10 


10 


(0,2,8) 


-0.093 


Z12 


12 


(0,0,12) 


-0.169 


71 A 


1 A 


(U,U,l/,/J 


n no*} 
-u.uyo 


Z15 


15 


(0,0,12,3) 


-0.037 


Z16 


16 


(0,0,12,4) 


+0.013 




FIG. 4: Distribution of We around Fe atoms, in supercooled Fe, FeB and FeZrB at T=800K 
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C. Bond Orientation Order Parameter We 

D 

Steinhardt, et al. 7] introduced the W\ parameters as a measure of the local orientational order in liquids and 
undercooled liquids. To calculate Wi, the orientations of bonds from an atom to its neighboring atoms are projected 
onto a basis of spherical harmonics. Rotationally invariant combinations of coefficients are then averaged over many 
atoms in an ensemble of configurations. The resulting measures of local orientational order can be used as order 
parameters to characterize the liquid structures. For an ideal icosahedral cluster, / = 6 is the minimum value of I for 
which Wq 7^ 0. Table [Til enumerates We values for different ideal clusters. The ideal icosahedral value of We is far 
from other clusters, making it a good icosahedral order indicator. 



Kasper-Polyhedrons 4, 38, 39] are polyhedrons which minimize the number of disclinations for a particular coor- 
dination number. The series Z8-Z16 in Table. |TTJ are such Kasper-Polyhedrons but with the added constraint that 
the surface atoms be triangulated with equilateral triangles. The icosahedron with a coordination of 12 is one such 
Kasper-Polyhedron with no disclinations. Adding disclinations to the icosahedron, one finds that each disclination 
increases the We value by the same amount irrespective of its sign. 

As before, we choose the cutoff distance to specify near neighbors as R cu t = 3.5A. For pure Fe at high tempera- 
tures, We resembles that of h igh temperature liquid Cu 1(, which in turn resembles the maximally random jammed 



configuration of hard spheres [24((Fig.[3]). On this basis we suggest that the MRJ configuration represents an idealized 
structure that is universal for strongly repulsive interactions. All pure metallic systems might approach this ideal 
structure at sufficiently high temperature. 

However, as temperature drops the We distribution shifts strongly to the left, with a pronounced increase in We < 
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-0.1 (Fig. [4]). This indicates a rather high concentration of nearly icosahedral clusters in supercooled liquid Fe. Even 
at very low percentages of supercooling of Cu (~ 3% l|) and Fe (~ 2%, not shown), Fe shows a clear enhancement 
in the negative Wq distribution as compared to Cu. 

We checked the W§ distribution of a non-magnetic simulation to see if it is strongly influenced by magnetism, and 
found a nearly identical result. In particular, we still found a strong enhancement of the nearly icosahedral clusters 
relative to liquid Cu or pure Fe at high temperature. 

D. Voronoi Analysis 



To explain the origin of this low We peak in pure Fe, we performed a Voronoi analysis 40j of the liquid before and 
during crystallization. A Voronoi polyhedron is described by indices (F3, F4, F$, ...) where fj denotes the number of 
faces with i edges. For example (0,0,12) denotes an icosahedron, while (0,0,12,2) denotes a 14-coordinated (Z14) atom, 
with 12 5-fold bonds and 2 6-fold bonds. The (0,0,12,2) is a characteristic TCP (tetragonal close-packed) structure 
of the Frank-Kasper type, with a —72° disclination line running through an otherwise perfect icosahedron. In a 
body-centered cubic crystal all atoms are of Voronoi type (0,6,0,8), which is an alternate 14-coordinated structure. 

Supercooled liquid Fe at T=800K (Table [TTT]) contains a high fraction of icosahedral atoms of type (0,0,12) and 
(0,0,12,2) . Those with Voronoi type (0,2,10,(0,1}) and (0,1,10,2) also have very negative We, so that they can be 
thought of as related to the icosahedron. Together, they explain the enhanced negative W§ distribution. Supercooled 
Cu (which shows weak icosahedral order) and high temperature Fe (T=1900K), in contrast contain about 1.2% of 
(0,0,12) and no (0,0,12,2). The MRJ configuration contains 1.2% (0,0,12) but no (0,0,12,2). 

Strikingly, the icosahedral clusters tend to join in pairs and strings of 3 atoms in length. An instantaneous quench of a 



TABLE III: Voronoi statistics (% of Fe atoms) for typical Fe-centered clusters in Fe, FeB and FeZrB. Edges smaller than 
0.33A and faces smaller than 0.3 A 2 have been treated as a single vertex. 



Voronoi 


Fe 


FeB 


FeZrB 


type 


supercooled 


defective 


supercooled 


supercooled 




liquid 


crystal 


liquid 


liquid 




800K 


relaxed 


800K 


relaxed 


800K 


800K 


(0,0,12) 


8.3 


11.0 


5.0 


0.0 


0.6 


1.4 


(0,0,12,2) 


2.1 


3.6 


2.0 


2.0 


0.5 


0.3 


(0,2,10,{0,1}) 


6.2 


6.8 


2.0 


0.0 


4.4 


7.1 


(0,1,10,2) 


3.6 


5.4 


3.0 


0.0 


2.1 


3.1 


(0,4,8,{0,1,2,3» 


7.9 


4.8 


10.0 


6.0 


8.6 


9.4 


(0,3,8,{0,1,2,3» 


6.8 


6.0 


12.0 


0.0 


4.3 


8.3 


(0,2,8,11,2,3,4}) 


3.0 


3.0 


3.0 


0.0 


0.0 


0.0 


(0,3,6,4) 


2.1 


2.2 


2.0 


6.0 


0.0 


0.0 


(0,5,4,4) 


1.5 


1.2 


3.0 


18.0 


0.0 


0.0 


(0,6,0,8) 


0.0 


0.4 


0.0 


56.0 


0.0 


0.0 



particular liquid structure at 800K which had a high fraction of icosahedral units, using conjugate gradient relaxation 
of the atomic coordinates and lattice parameters, shows a clear enhancement of the icosahedral and other closely 
related units. The strings of icosahedral units found in the supercooled liquid connected to form networks. A quench 
starting from a different instantaneous supercooled liquid structure at 800K containing fewer icosahedral structures 
resulted in rapid crystallization. This indicates that the presence of icosahedrons may inhibit crystallization. Similar 
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FIG. 5: Fragment of BCC crystal illustrating a central bond in the [1,1,1] (vertical) direction, surrounded by two equilateral 
triangles rotated 60° from each other and displaced on either side of the midplane. Balls and sticks follow ZomeTool convention: 
Balls are icosahedral, yellow sticks are 3-fold and blue sticks are 2-fold bonds. 

quenches, starting from instantaneous liquid structures at higher temperatures, also resulted in relaxed structures 



with some of them showing a high negative We distribution comparable to the distribution at 800K. The quenches 
that partially crystallized showed a high fraction of BCC (0,6,0,8) 's (~ 40%) and was always accompanied by a high 
fraction of (0,5,4,4)'s (~ 25%), which are otherwise absent or very low in the liquid. 

Under TMD the supercooled liquid at T=800K eventually crystallizes, with the W§ distribution peaked strongly 
around zero consistent with the value for ordinary crystalline clusters, but retaining a subset of atoms with nearly 
icosahedral Wq<-0.14:. Not surprisingly, these were precisely the atoms that had Voronoi type (0,0,12) prior to 
crystallization. The Z14 (0,0,12,2) atoms (with We ~ -0.093) were mutual near-neighbors, linked along their 6-fold 
(—72°) bonds, and also were neighbors of the nearly icosahedral atoms. 

We quenched this sample by conjugate gradient relaxation of atomic coordinates and lattice parameters. After 
relaxation we found 56 atoms had BCC Voronoi type (0,6,0,8). Six icosahedral atoms became 12-coordinated (0,4,8) 
structures surrounding the bond connecting the two (0,0,12,2) disclinated icosahedral atoms, which retained their 
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type. The remaining atoms served to link the cluster of icosahedron-related atoms to the surrounding defect-free 
BCC crystal. 

The icosahedron-related structure thus forms a point defect in an otherwise perfect BCC crystal. A simple way to 
create this defect is to take two consecutive triangles surrounding a near-neighbor bond along the BCC [1,1,1] direction 
(see Fig. [5]), displace them into the perpendicular bisecting plane and rotate by 30°. The atoms along the [1,1,1] bond 
are now connected to each other by a six-fold bond (—72° disclination) and are connected to the six displaced atoms 
by five-fold bonds. This scheme to transform a BCC rhombic dodecahedron into a Frank-Kasper Z14 polyhedron was 



described in Ref 4l|, |42j in an attempt to explain certain diffraction anamolies in iron and vanadium-based alloys 



under ion irradiation. 

Manually removing the defect, by reversing the above procedure then relaxing, yields a perfect BCC crystal with 
all atoms in a (0,6,0,8) Voronoi environment. We also embedded the point defect in an otherwise perfect BCC crystal 
with 128 atoms in a cubic box. Relaxation showed that the defect was stable. The energy of the defect was 6.18eV 
and the fractional volume increase was 0.013. 

Since tungsten, like Iron, crystallizes in BCC, we performed a separate first-principles simulation of liquid Tungsten. 
After supercooling by 13% to T=3200K, we found the We histogram resembled that of Fe at 1600K (which is 15% 
supercooled). A Voronoi analysis revealed a similar percentage of (0,0,12)'s and (0,0,12,2)'s in both these elemental 
BCC-forming metals. Hence we believe the relationship between BCC and icosahedral structures may be linked to 
our observed high fraction of icosahedra in liquid Fe as compared to Cu. This may also explain the reason why we 
are able to supercool BCC Fe more deeply than FCC Cu in our simulations. 

We also included an icosahedron point defect inside an otherwise perfect FCC crystal of Cu, with 256 atoms in a 
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FIG. 6: Honeycutt- Andersen analysis for pure Fe shows a clear increase in five-fold bonds with supercooling. 



cubic box. Relaxation showed that the defect was stable. The energy cost of the defect was 4.69eV and the fractional 



volume increase was 0.012. Even though this FCC defect cost less energy than the point defect in BCC, and also needs 
little rearrangement of atoms as reflected by the slightly lower fractional density change, we do not see a significant 
icosahedral order in liquid Cu, when compared to liquid Fe. This could be because the number of icosahedrons that 
a single Z14 disclination can stabilize (up to six) is greater than one. So even an equal number of the two different 
defects in BCC and FCC would result in more icosahedral order in BCC than in FCC. 



Honeycutt and Andersen [6( introduced a useful assessment of local structure surrounding interatomic bonds. We 
employ a simplified form of their analysis, counting the number of common neighbors shared by a pair of near-neighbor 
atoms. This identifies the number of atoms surrounding the near-neighbor bond and usually equals the number of 



edge-sharing tetrahedra whose common edge is the near-neighbor bond. We assign a set of three indices to each 



E. Honeycutt and Andersen analysis 



fl 
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bond. The first index is 1 if the root pair is bonded (separation less than or equal to R cu t)- The second index is 
the number of near-neighbor atoms common to the root pair, and the third index gives the number of near-neighbor 
bonds between these common neighbors. We take the same value of R cu t=3.5A as mentioned before. Note that the 
Honcycutt and Andersen fractions depend sensitively on R cu t, making precise quantitative comparisons with other 
prior studies difficult. 

In general, 142 's are characteristic of close packed structures (FCC and HCP) and 143 's are characteristic of 



distorted icosahcdra 



43j. They can also be considered as +72° disclinations [3j,[4], [5|. Likewise, 15's are characteristic 



of icosahedra, with 155's characterizing perfect icosahedra while 154's and 143's characterize distorted ones. 16's 
indicate —72° disclinations. 166's and 144's are also characteristic of BCC. 

A Honeycutt and Andersen analysis for pure Fe, with an R cut —3.5A (Fig. [6]), showed that with supercooling the 
fraction of 15 bonds rises from 0.46 at T=1900K to 0.59 at T=800K in the liquid before crystallization. The fraction 
of 155's (characteristic of perfect icosahedra) was always larger than 154's (characteristic of distorted icosahedra), and 
seemed to be steeply increasing with supercooling as opposed to 154's, which were relatively flat. The fraction of 14 
bonds drop from 0.32 to 0.30, with the icosahedral 143's being always higher than the cubic 142's. The 144's, which 
are characteristic of BCC remain nearly flat, even though the 166's show a slight increase. The ease of embedding a 
Z14 disclination in BCC Fe (see section [ill Dl) might explain the slight increase in the 166's. 



19 




FIG. 7: Partial radial distribution functions in FeZrB at T=800K. 
IV. FE-B AND FE-ZR-B 

A. Radial Distribution Function g(r) 

Fig. shows the pair correlation functions of supercooled FeZrB at T=800K. From the heights of the first peaks, 
we see that the strongest bonds form between the metalloid (B) and the metal (Fe or Zr). The relative bond lengths 
reveal, as expected, that B behaves as a small atom, Fe is medium sized and Zr is large. Comparing simulations of 
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TABLE IV: Average coordination number in FeB. 





FeFe 


FeB 


BFe 


BB 


Rcut 


3.4 


3.0 


3.0 


2.3 


N af} 


12.0 


2.2 


8.9 


0.4 



TABLE V: Average coordination number in FeZrB. 





FeFe 


FeZr 


FeB 


ZrFe 


ZrZr 


ZrB 


BFe 


BZr 


BB 


Rent 


3.4 


3.8 


3.0 


3.8 


4.0 


3.3 


3.0 


3.3 


2.4 




9.5 


2.0 


2.0 


14.0 


0.9 


2.8 


6.8 


1.4 


0.4 



FcZrB to pure Fe, we find that gFeFe is reduced by alloying with B or ZrB, because Fe prefers to associate with B or 
Zr rather than with Fe. 



To calculate the coordination number, define N a p as the average number of atoms of type /3 around an atom of type 
a, We set R cu t at the first minima of the partial radial distribution functions (Fig. [7]). We list the partial coordination 
numbers of FeB in Table. IIVI and FeZrB in Table. [V] (averaged over all the temperatures since the temperature 
dependence is very weak and non-monotonic). The average value of Np e Fe decreases with alloying, due to decrease in 
the concentration of Fe and the preference to bind with B and Zr. Zr, being a large atom, has a larger coordination 
number. Also note that B and Zr favor each other more than themselves. We do find some B-B pairs in the liquid 



state. 
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FIG. 8: Simulated partial structure factors of FeB at T=800K 



B. Liquid Structure Factor S(q) 



Fig.[S]shows the Faber-Ziman 44( partial structure factors of FeB at the lowest simulated temperature of T=800K, 



defined as, 



S af3 (q) = 1 + 47rp / [g afJ (r) - 1]^* 



qr 



(4) 



The positions of the first and second peaks in partial SFeFe (<?) is in very good agreement with the experimental results 



for amorphous FeB 45j. Also, at the position of the splitting of the second peak in SVeFe(<7), as observed in the 



experiments, we observe a slight shoulder. Similarly, the positions of the different peaks in SFeB(q) and Sbb(q) are 



in agreement with the experiments. 



For FeB, the q — > limit of SFeFe(q) is comparable to the experimental value. The q — > limit in the other two 
partials differ from the experiment, SsB(q) more seriously than SFeB{q)- We think that this discrepancy in the long 
wavelength regime is probably due to the very low density of B in our system leading to poor statistics. Nevertheless, 



the excellent agreement in the positions of the different peaks in the partial structure factors shows that we have 
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FIG. 9: Distribution of We for different chemical species in supercooled FeZrB at T=800K 



reasonably good representative structures of FeB at T=800K. Partial structure factors or partial pair distribution 



functions are not available experimentally to compare with FeZrB simulations. 



C. Bond Orientation Order Parameter We 



To define the We distribution in an alloy, we concentrate on central atoms of some particular species (e.g. Fe) but 
consider the neighboring atoms of all species. We chose the near-neighbor cutoff distances as before. Compare the 
Fe-based We distributions at their supercooled temperatures in Fig. [4] The origin of high negative We values for pure 
Fe was previously explained in section IIIIDl Replacing a few medium sized Fe atoms with smaller B atoms causes 
negative disclination lines to concentrate on Fe, leading to a drop in the ideal icosahedral clustering on these Fe atoms 
and strongly reducing the extreme negative values of We in FeB. On the contrary, inclusion of Zr in FeZrB causes 
negative disclinations to attach to Zr, easing frustration, leading to more Fe centered clusters with icosahedral ordering 



and increasing the negative region of We, as compared to FeB. Inclusion of the big Zr atom enhances icosahedral 
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order on the Fe atoms. 

Fig. [5] shows the W§ distributions for FeZrB with centers at Fe, Zr and B. The histogram with the center at Zr 
is almost symmetric about the value of zero. This suggests that the local environment about Zr atoms is nearly 
spherical, as is expected given its large size. The B centered Wq histogram is also asymmetric towards negative Wq 
values due to Kasper-Polyhedrons and slightly distorted versions of them (see Table. HI1 and Section fill Cp . 

D. Voronoi Analysis 

A Voronoi analysis was performed for FeB in the supercooled liquid at T=800K (Table [TTT|) . The Fe environments 
were mostly (0,4, 8,x) or (0,3,8,x) types, where x={0,l,2,3,4}, with the higher coordination polyhedron being more 
favored. This was followed by the very negative Wq valued (0,0,12), (0,2,10,x)'s and (0,l,10,2)'s Voronoi types 
occurring at lower frequency than in pure Fe. Boron mainly had environments of type (0,3,6) 's (~ 15% of B atoms) 
(Kasper polyhedron for Z=9 containing a +72° disclination) and (0,5,4) 's (~ 20%). These types are typical of the 
tri-capped trigonal prism (TTP) and the mono-capped square archimedean prism (a slightly distorted variant of 
TTP) respectively. The TTP is found in the crystal structures of FesB with Pearson symbols oP16 and tI32. The 
distorted (0,5,4) version is found in the structure Fe23B6 of Pearson type cF116. These structures have been identified 



as the leading competitors for B-Fe glass [46[. Boron also took environments of the Kasper polyhedron (0,4,4) (~ 
3%) corresponding to Z=8 and (0,2,8) (~ 10%) corresponding to Z=10. The association of B with +72° disclinations 
explains how it increases the frustration of icosahedral order on the Fe atoms. Clearly the improved glass-formability 
of FeB compared with elemental Fe cannot be due to icosahedral order. Rather, it is presumed to be caused by the 
deep eutectic at Fes3Bi 7 . 
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A Voronoi analysis of supercooled FeZrB at T=800K shows a clear increase in the very negative Wq polyhedrons, 
and also a decrease in the number of Z14 (0,0,12,2) types on Fe atoms indicating a decrease in frustration in the 
ternary as compared to the binary. Environments around B atoms were roughly similar in the binary and the ternary, 
with a slight increase in the lower coordinated (Z=8) Kasper polyhedrons at the cost of higher coordinated (Z=10) 
ones. Zirconium took a variety of polyhedra, with an average coordination of 17.6 and a minimal coordination of 15, 
owing to its large size compared to the other constituents. 



E. Honeycutt and Andersen analysis 



We made a HA analysis of the ternary glassy alloy, by looking at root pairs of chemical species a and (3, choosing 
R cu t i n the manner of the Wq analysis. The frequency is normalized to sum to one for each species pair a/3. Table I VII 
lists the fraction of different Ix pairs in supercooled FeZrB. 

Among the lx pairs with Fe as one of the root pairs, the 15's are most abundant at all temperatures. The 15's 
are mainly comprised of 155's and the 154's, with the 155's being always higher than the 154's. The percentage of 
15's is similar for FeZrB, FeB as well as pure Fe. Note that the 15's are largest for the Fe-Fe pairs. They also show 
a steady enhancement with supercooling, unlike the 14's which decrease with supercooling. Among the 14's, the 
icosahedron-related 143's for all root pairs are always higher than the close-packed 142's, and remain fairly constant 
with supercooling. The 14's are maximal for FeB root pairs and minimal for FeZr. 

Supercooled pure Fe has a high percentage of "16" 's (13%) compared to pure Cu (7%). The high number of "16" 's 
in Fe is related to the occurance of Z14 (0,0,12,2) environments in which the 6-fold bond carries a —72° disclination, 
rather than the BCC 14 atom arrangement which would also show a high degree of four-fold and six-fold bonds 



25 



TABLE VI: HA analysis for supercooled FeZrB at T=800K 



Root Pair 




Fe-Fe 


Fe-Zr 


Fe-B 


Zr-Zr 


Zr-B 


B-B 


14 pairs 


0.30 


0.19 


0.54 


0.11 


0.55 


0.75 


142 pairs 


0.08 


0.06 


0.03 


0.06 


0.03 


0.02 


143 pairs 


0.17 


0.11 


0.35 


0.04 


0.33 


0.33 


144 pairs 


0.05 


0.02 


0.17 


0.01 


0.19 


0.41 


15 pairs 


0.56 


0.53 


0.39 


0.48 


0.41 


0.18 


154 pairs 


0.22 


0.27 


0.03 


0.36 


0.04 


0.0 


155 pairs 


0.33 


0.22 


0.37 


0.08 


0.36 


0.17 


16 pairs 


0.10 


0.24 


0.01 


0.36 


0.01 


0.01 


166 pairs 


0.09 


0.20 


0.0 


0.26 


0.00 


0.00 



( section IlII D[) . Adding B to Fe, shows an increase of "16" 's for Fe-Fe pairs to 21% in FeB (not shown). Adding 
large Zr atoms to FeB decreases the occurance of 16's on the Fe-Fe pairs, putting them on Fe-Zr pairs and easing the 
frustration of Fe centers. This causes the geometry about Fe centers to be more icosahedral, and hence the shift of 



the We towards negative values. 

Among the root pairs not containing Fe, note that the B-B pair has the maximum 14's (+72° disclinations) , 
especially 144's, while the Zr-Zr pair has the maximum 16's (—72° disclinations), especially the 166's, emphasizing 
the role of size in controlling the frustration in alloys. 
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V. CONCLUSION 

This study quantifies icosahedral and polytetrahedral order in supercooled liquid metals and alloys. This is the 
first such analysis of glass-forming Fe compounds using configurations from first-principles simulations. While the 
structural properties of Fe and Cu strongly resemble each other at high temperature, and indeed are close to a 



maximally random jammed structure l|, |24| |. their behavior evolves substantially, and in different manners, as the 
liquid is supercooled. Proper modeling of atomic interactions is essential to capture the differing behavior of each 
element, and use of a first-principles simulation is the most reliable means of achieving this. 

For pure elements we find the degree of local icosahedral order in the supercooled liquid depends on the low 
temperature crystal structure, with BCC metals such as Fe and W accommodating icosahedra more readily than the 
FCC element Cu. Alloying with large or small atoms can further influence the degree of icosahedral order, with small 
atoms (e.g. B in Fe) aggravating the frustration by introducing positive disclination line defects, while large atoms 
(e.g. Zr in Fe) naturally stabilize negative disclination line defects, relieving frustration on the medium-sized Fe atoms. 
While the enhanced glass-forming ability of FeB compared to Fe cannot be related to icosahedral order, we suggest 
that the formation of icosahedral order and disclination line network, together with the slow dynamics of chemical 



ordering in a complex alloy and the destabilization of competing crystal phases 46j , enhances the glass-forming ability 
of FeZrB compared with FeB. 
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